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I present an introduction to tlie NRQCD factorization mctiiod for calculating an- 
nihilation rates and inclusive production rates for heavy quarkonium. Using this 
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X \ • method, annihilation decay rates and sufficiently inclusive cross sections are fac- 
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tored into long-distance NRQCD matrix elements and perturbative short-distance 
coefficients, f derive the velocity-scaling rules that are used to estimate the magni- 
ff^ , tudes of the nonperturbative NRQCD matrix elements and I describe perturbative 

matching methods for calculating their short-distance coefficients. Some simple 
applications for which NRQCD factorization methods have dramatic implications 
are discussed. 
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p^ ' 1 Introduction 

(N : 

f ^ I The NRQCD factorization approach is a systematic framework for analyzing 

t~^ . annihilation decay xates and sufficiently inclusive production cross sections for 

0^ I heavy quarkoniumtl This method separates the effects of short distances that 

are comparable to or smaller than the inverse of the heavy quark mass from 
Qh| the effects of longer distance scales. Short-distance effects are calculated using 

perturbative QCD and long-distance effects are described by matrix elements 
fl_j ' in an effective field theory called nonrelativistic QCD (NRQCD). This paper is 

<~i . an introduction to the NRQCD factorization approach. In Section pi the effec- 

tive lagrangian for NRQCD is constructed, the velocity-scaling rules that are 
used to estimate the magnitude of NRQCD matrix elements are derived, and 
perturbative matching methods for calculating short-distance coefficients are 
described. The application of NRQCD factorization methods to annihilation 
rates and to inclusive production rates are discussed in Sections g and ]^ 

2 Heavy Quarkonium and NRQCD 

2.1 Scales in Quarkonium Physics 

Heavy quarkonium is a meson containing a heavy quark and its antiquark. 
There are a number of different energy scales that play an important role in 
the quarkonium physics. While the many scales make the physics complex, 
they also make it interesting. The scales of quarkonium include the mass 
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Table 1: Quarkonium energy scales 
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bb 


tt 


M 

Mv 

Mv2 


1.5 GeV 
0.9 GeV 
0.5 GeV 


4.7 GeV 
1.5 GeV 
0.5 GeV 


180 GeV 
16 GeV 
1.5 GeV 



M of the heavy quark, its typical momcntuni Mv, and its typical kinetic 
energy Mw^. The quark mass M sets the total energy scale for annihilation 
decays and the scale of the kinematic threshold for onium production. The 
inverse of the typical momentum Mv is the length scale for the size of the 
onium state. The typical kinetic energy Mv'^ is the scale of the splittings 
between radial excitations and between orbital-angular-momentum excitations 
in the onium spectrum. For both charmonium and bottomonium, the splittings 
between the two lowest '^S'l states ( J/V' and V'' for charmonium, T and T' for 
bottomonium) are approximately 600 MeV. The splittings between the lowest 
^Si and ^Pj states {J/ip and Xcj for charmonium, T and Xbj for bottomonium) 
are approximately 400 MeV. We take 500 MeV as an estimate for the scale 
Mv"^. In Table 0, we give our estimates for the scales M, Mv, and Mv^ 
for charmonium and bottomonium. The estimate for M is half the mass of 
the lowest energy level of quarkonium. The estimate for the scale Mv is the 
geometric mean of the estimates for M and Mv"^. From Table ^ we see that v'^ 
is approximately 1/3 for charmonium and 1/9 for bottomonium. These values 
are small enough to justify a theoretical approach based on an expansion in 
powers of v. 

Another important energy scale in quarkonium physics is Aqcd , the scale 
of nonperturbative effects involving gluons and light quarks. The potential 
energy V{R) between a quark and antiquark separated by a distance R varies 
from Coulombic at sufficiently short distances, 
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to linear at sufHciently long distances: 
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where empirically k ~ 
proportional to Aqcd- 



450 MeV. Since k is independent of M, it must be 



To clarify the relation between the scale Aqcd and the other scales in 
quarkonium physics, we consider two limiting cases. The first case is a heavy 
quark whose mass M is large enough that the onium wavefunction is dominated 
by the Coulombic term in the potential. The size of a bound state is determined 
by a balance between the kinetic energy and the potential energy, so we must 
have 
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The size of the bound state is comparable to the inverse of the typical momen- 
tum Mv of the heavy quark in the bound state. Setting R ^ l/{Mv) in (0), 
we obtain 

V ~ as{Mv). (4) 

This equation can be solved self-consistently for u as a function of M. If M is 
sufficiently large, the resulting value of Mv'^ is much greater than Aqcd- For 
the top quark with mass Af = 180 GeV, we obtain v ~ 0.18. An alternative 
estimate for v, which is closer in spirit to the estimates used for charmonium 
and bottomonium in Table ||, is obtained by taking Mv"^ to be the splitting 
between the ground state and the first excited state in the potential (|l|) . This 
splitting is approximately 1.5 GeV, which corresponds to u ~ 0.09. This value 
for Mv'^ and the corresponding estimate for the scale Mv are included in 
Table H. These would be the appropriate scales for toponium if the top quark 
were stable enough to form bound states. 

The second limiting case is a heavy quark whose mass M is small enough 
that the wavefunction is dominated by the linear term (0) in the potential, but 
still large enough that the onium is nonrelativistic. The balance between the 
kinetic and potential energies requires 

Mv^ - k'^R. (5) 

Setting R ~ l/{Mv), we find that k ~ Afw'^/^. Identifying k with the scale 
A-QCD, we find that the typical velocity of the heavy quark is such that this 
scale is intermediate between the scales Mv and Mv"^ . If the Coulombic and 
linear regions of the potential were both equally important, then the relations 
(^ and (|^) would be satisfied simultaneously. The fact that these relations are 
compatible suggests that the scaling relation (Q) might be applicable even for 
rather low values of M . 

It is a remarkable coincidence of quarkonium physics that the scales Mv'^ 
are almost identical for charmonium and bottomonium. The fact that this scale 
is so insensitive to the value of M suggests that Aqcd should be identified with 
the scale Mv"^ for quark masses in this range. This identification is supported 
by the numerical value of the scale Mv'^ given in Table |[ 



Table 2: Value of the QCD coupling constant at the characteristic momentum scales for 
heavy quarlionium 





cc 


bb 


ti 


as{M) 
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0.35 
0.52 
> 1 


0.22 
0.35 
> 1 


0.11 
0.16 
0.35 



If any of the scales M, Mv, and Mv"^ is large enough compared to Aqcd, 
then the effects of that scale can be calculated using QCD perturbation theory. 
The values of the running coupling constant of QCD at the scales M, Mv, and 
Mv"^ are given in Table y for charmonium, bottomonium, and toponium. All 
three scales are perturbative in the case of toponium. For charmonium and 
bottomonium, as{M) is small enough to justify perturbation theory at the 
scale M. The coupling constant as{Mv) is also small enough that perturbation 
theory seems reasonable. However, potential model calculations imply that the 
wavefunctions of charmonium and bottomonium have significant support in the 
linear region of the potential, and this suggests that nonperturbative effects 
may be significant at the scale Mv. As for the scale Mw^, it is hopelessly 
nonperturbative for charmonium and bottomonium. 

In addition to the scales M, Mv, Mv"^, and Agcu, there are also kine- 
matic energy scales that can play an important role in quarkonium physics. 
For example, in the production of quarkonium, the total center-of-mass energy 
Y^ and the transverse momentum pT of the onium can be important. Unrav- 
elling the effects of the various energy scales is essential in order to understand 
quarkonium physics. This is particularly important for charmonium and bot- 
tomonium, because the coupling constant as runs rather dramatically with the 
momentum at scales of order M and smaller. 

The NRQCD factorization approach is based on separating short-distance 
effects involving momenta of order M or larger from those effects that involve 
the smaller momentum scales Mv, Mv"^, and Aqcd- The scale M is assumed 
to be perturbative, so that short-distance effects can be calculated using per- 
turbation expansions in as{M). No assumption is made about the validity of 
perturbation theory at the scale Mv. Instead, we exploit the fact that in a 
nonrelativistic bound state, the typical velocity v provides a small expansion 
parameter. 



2.2 Integrating out Relativistic Effects 

One way to separate the effects of the momentum scale M from the lower 
momentum scales in a field theory is to integrate out all modes with momenta 
greater than some cutoff A that is much less than M. The result of this 
renormalization group transformation is a Wilsonian effective field theory that 
describes the modes with momenta smaller than A. All effects of the scale 
M are encoded in the parameters of the effective field theory. In our case, 
the original field theory is QCD with a heavy quark. It is described by the 
lagrangian 

£qcd = Cuei,t+^{i7^D^-mQ)^, (6) 

where >Ciight is the lagrangian that describe gluons and light quarks. The mass 
parameter mq of the heavy quark can be identified with M . It is implicit 
in the Lorentz-invariant lagrangian m) that the cutoff is much larger than 
M. Integrating out the momentum scale M is equivalent to lowering the 
cutoff A to a value lower than M. We will argue that if A is in the range 
Mv <C A ^ M, the resulting effective field theory can still be described by a 
local lagrangian. In other words, the effects of modes with momenta of order 
M can be reproduced by local interactions among the lower momentum modes. 

Suppose the onium is in a virtual state that includes a quark with rela- 
tivistic momentum of order M . Then that state is off its energy-shell by an 
amount of order M, which is much larger than the scale AIv'^ of the splittings 
between onium energy levels. By the uncertainty principle, the lifetime of that 
highly virtual state is less than or of order 1/M. In that short time, fields 
can propagate only over distances of order 1/M that are pointlike on the scale 
l/{Mv) of onium structure. Thus the effects of virtual states that are excluded 
by a momentum cutoff in the range AIv <^ A <ti M can be reproduced by local 
interactions involving low- momentum modes. 

The above argument applies equally well to virtual states that contain a 
light parton with momentum of order M in addition to the QQ pair. However, 
it does not apply to virtual states obtained by the annihilation of the QQ pair. 
Such states can contain light partons with momenta of order M without being 
far off the energy shell. A momentum cutoff satisfying A ^ M excludes these 
states, but their effects cannot be reproduced in detail by local interactions. 
For example, the annihilation decay of the onium produces light hadrons, some 
of which must have momenta of order M and therefore parton constituents with 
momenta of order M. With such modes excluded by the cutoff, we cannot hope 
to describe the annihilation decays accurately. Nevertheless, as we shall see 
in Section 0, the total annihilation width of an onium state can be described 
accurately. Thus the momentum cutoff A can be extended to light partons at 



the expense of a restriction on the physical observables that can be described 
within the effective theory. 

With a momentum cutoff A that excludes relativistic Q and Q states, it is 
convenient to describe the heavy quark and antiquark by separate 2-component 
Pauli fields -0 and x, rather than by a single 4-component Dirac field \1/. ff we 
simply substitute "^ — y^ into the lagrangian (|g), we obtain off-diagonal terms 

that couple ^^ to x and x^ fo V'l allowing the creation and annihilation of QQ 
pairs. We will argue that terms that change the numbers of heavy quarks and 
antiquarks can be removed from the lagrangian and compensated by terms 
that conserve the numbers of Q's and Q's. This is obvious for terms that 
create QQ pairs, because a virtual state containing an additional QQ must be 
off its energy shell by an amount of order M . As mentioned above, this is not 
completely true for terms that allow the QQ pair in the onium to annihilate. 
The effects of states consisting of gluons and light qq pairs that are produced by 
QQ annihilation cannot be reproduced in detail by local interactions. However, 
the effects of these states on sufficiently inclusive observables can be described 
accurately. Thus, with this restriction on physical observables, terms in the 
effective lagrangian that change the numbers of Q's and Q's can be eliminated 
from the effective lagrangian. 

For the Dirac term of the lagrangian (^, the decoupling of the fields ip 
and X can be accomplished by a unitary transformation called the Foldy- 
Wouthuysen-Tani transformation. For the case of a background gauge field, 
it is straightforward to construct the transformation that diagonalizes the la- 
grangian to any desired order in the heavy quark velocity. The simplest form 
of this transformation in the Dirac representation is 

* -^ exp(-i7-D/2mQ)^. (7) 

After this transformation, the heavy-quark term in the lagrangian (0) can be 
approximated by 

0\ t / -mg + iDo + T>y2mQ \ f^\ , . 

Xj \ mQ+iDo--Dy2mQ jyxj' ^' 

If we take D to scale like Mv, the corrections to the entries in the matrix 
scaling like Mv"^. 

For the effective lagrangian with momentum cutoff A, the elimination of 
terms that change the numbers of Q's and Q's is more complicated than simply 
applying a unitary transformation. In addition to the terms that are quadratic 
in ip and Xj the effective lagrangian also includes terms that are quartic and 

6 



higher in the heavy quark fields. A further comphcation is that gluon interac- 
tions modify the coefficients of the terms produced by the Foldy-Wouthuysen- 
Tani transformation. Nevertheless, by the general arguments presented above, 
one can describe the low-energy QQ sector of QCD by an effective field theory 
in which the numbers of heavy quarks and antiquarks are strictly conserved. 

2.3 Effective Field Theory 

One could in principle construct a nonrelativistic effective lagrangian that de- 
scribes the low-energy QQ sector of QCD by starting with the lagrangian (ra) 



and carrying out the sequence of two transformations described in Section 2.2 
The first is a renormalization group transformation that removes modes with 
momenta greater than A. The second is a transformation that removes interac- 
tions that change the numbers of heavy quarks and antiquarks. Both of these 
steps would be extremely complicated to carry out in practice. Fortunately, 
there is an alternative to the explicit construction of the effective lagrangian 
and that is to use the strategy of "effective field theory" 13 In this approach, the 
construction of the effective lagrangian proceeds through the following steps: 

1. Identify the fields that are required to describe the low-energy excitations 
of the theory. 

2. Identify the symmetries that one could maintain in the effective theory 
by using a suitable cutoff and making appropriate field redefinitions. 

3. Specify the accuracy to which low energy observables in the original 
theory should be reproduced by the effective theory. 

4. Write down the most general effective lagrangian that is consistent with 
the symmetries, including all terms that are required to reproduce the 
physics to the specified level of accuracy. 

5. Determine the coefficients of those terms by matching low-energy observ- 
ables of the effective theory with those of the full theory. 

The effective field theory that is obtained by applying the above strategii. 
to the low-energy QQ sector of QCD is called nonrelativistic QCD (NRQCD)El 
The fields that are required to describe the low energy degrees of freedom are 
the heavy quark and antiquark fields V' and x, the S'[/(3) gauge fields A^, 
and the Dirac fields for the light quarks. The symmetries of NRQCD are the 
following: 



• 



• 



• 



• 



SU{3) gauge symmetry. This local symmetry requires that the gluon 
fields enter into the effective lagrangian only through the gauge- covariant 
derivatives Dq and D and the QCD field strengths E and B. 

rotational symmetry. A nonrelativistic description of the heavy quark 
necessarily breaks the Lorcntz symmetry of QCD down to its rotational 
subgroup. 

charge conjugation and parity. These discrete symmetries of QCD are 
also symmetries of the effective theory. The charge conjugation transfor- 
mations of the heavy quark and antiquark fields are 

V' ^ »(X^^2)*, X ^ -«(V'V2)*. (9) 

The parity transformations are 

V^(i,r) ^ ^(i,-r), x(i,r) -> -x{t,-r). (10) 

heavy-quark phase symmetry. This symmetry guarantees the separate 
conservation of the number of heavy quarks and antiquarks. Its action 
on the fields is 

V- - e^^V-, X ^ e^'^X- (11) 

Having identified the symmetries of NRQCD, we can write down the most 
general effective lagrangian that is consistent with these symmetries. It has 
the form 

-Cnrqcd = flight + ^^ (iDo + ^) V + X^ (iDo ~ ^) x + S£, (12) 

where £iight is the usual lagrangian that describes gluons and light quarks. The 
desired level of accuracy is specified by the order in v with which the onium 
energy levels must be reproduced by the effective theory. The heavy quark 
terms that are shown explicitly in (n3) are those that are required to calculate 
the energy levels up to errors of order Mw*. The term 5C in ( [l^ ) includes the 
correction terms that must be added to decrease the errors to order Mv^ or 
smaller. 

The invariance of physical quantities under field redefinitions can be ex- 
ploited to eliminate some terms in the NRQCD lagrangian. Gauge invariance 
requires that the gluon fields in 5C appear only through the covariant deriva- 
tives Dq and D. However, a redefinition of the field tp can be used to eliminate 
terms in which Dq acts on i/;, and similarly for x- Thus field redefinitions 



can be used to eliminate all occurences of Dq in SC except in the combination 
[Dq, D] = zgE. Because of these field redefinitions, NRQCD will not reproduce 
the low-energy behavior of the Green's functions of QCD. It will only agree 
with full QCD for on-shell physical quantities. 

The minimal form of NRQCD is obtained by setting SC = in ([l2|). It 
contains two parameters, the heavy-quark mass parameter M and the gauge 
coupling constant g. These parameters can be tuned as functions of the QCD 
coupling constant as , the heavy-quark mass parameter toq , and the ultraviolet 
cutoff A of NRQCD so that the splittings between the onium energy levels are 
reproduced up to errors of order Mv*. Since the energy splittings between ra- 
dial excitations (such as J/ip and ip') and between orbital- angular- momentum 
excitations (such as J/ip and Xcj) scale like Mv^, these are reproduced up 
to errors of relative order v^ . Spin splittings in heavy quarkonium, such as 
the splitting between the lowest ^5*0 and ^Si states {jjc and J/-;/; for charmo- 
nium), scale like Mw^. These splittings vanish in minimal NRQCD due to the 
following symmetry: 

• heavy-quark spin symmetry. Under this symmetry, the two spin compo- 
nents of the heavy quark and the two spin components of the antiquark 
are mixed by independent unitary transformations: 

V' ^ C/V, X -^ Vx. (13) 

where U and V are SU(2) matrices. This is only an approximate sym- 
metry of the complete NRQCD lagrangian, holding up to corrections of 



If we wish to reduce the errors in the quarkonium energy levels to smaller 
than order Mw^, it is necessary to add additional terms 5C to the lagrangian 
in dl2). Using the velocity-scaling rules that are discussed in Section 



2.4 



it 

can be shown that the terms that are necessary and sufficient to reduce the 
errors to order Mv^ are 






^^ " ^V^t(D2)V + ^V't(D.gE-.gE.D)V' 



, g^fV-^l^D X .gE - .gE x zD) • a^ + ^^^(^B • a)^ 
+ charge conjugate terms, (14) 

where ci, C2, C3, and C4 are dimensionless coefficients. We will refer to the 
terms in ( |l4[ ) as the v^ -improvement terms in the NRQCD lagrangian. The 
two terms in (^) that contain the Pauli matrix a break the spin symmetry 



of minimal NRQCD. They give spin splittings that scale like Mv* and are 
accurate up to errors of relative order v'^ . Splittings between radial excitations 
and splittings between orbital-angular-momentum excitations are reproduced 
up to errors of relative order u^. 

2.4 Velocity-scaling rules 

The relative importance of the terms in the NRQCD lagrangian can be de- 
duced from the self-consistency of the quantum field equations.. for minimal 
NRQCD and from the basic qualitative features of quarkoniumB The results 
of this analysis are summarized by the velocity-scaling rules in Table 0. The 
magnitude of a matrix element of a local gauge-invariant operator between 
quarkonium states can be estimated by multiplying the appropriate factors 
from Table ^. The scaling with M follows simply from dimensional analysis, 
so the nontrivial content of Table H is the scaling with v. 

The first few lines in Table ^ can be derived very easily. The expectation 
value of the number operator }<Px tp''tp i^i & quarkonium state \H) is very 
close to 1: 

{H\ [ d^xtP^'tP \H) w 1. (15) 



We have normalized the quarkonium state so that {H\H) = 1. From the 
fact that a quarkonium state can be localized to within a region 1/ {Mv)'^ , we 
conclude that ip must scale like (M?;)^/^. The expectation value of the kinetic 
energy term in the NRQCD hamiltonian scales like Mw^ : 

{H\ f cPx'ilj^D^/2M)i;\H) - Mv^ . (16) 

This implies that D must scale like Mv. The fact that Dq scales like Mv'^ 
when acting on ip then follows immediately from the field equation for ip: 

^Do-^)^^0. (17) 

The estimates for gAo and gA in Table are specific to Coulomb gauge, 
which is defined by V • A = 0. As shown below, the field equations in this 
gauge indicate that the effects of the vector potential A are suppressed by a 
factor of V relative to the scalar potential Aq . The dominant terms in the field 
equations for ip and Aq are therefore 

ido-gA, + ^)^ « 0, (18) 

VVo+ffW « 0. (19) 
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Table 3: Estimates of the magnitudes of NRQCD operators for matrix elements between 
heavy-quarkonium states. 



Operator 


Estimate 


V; 


(Mw)3/2 


X 


(Mw)3/2 


Dq (acting on ?/> or x) 


Mv^ 


D 


Mv 


gE 


M^v^ 


gB 


M'^v^ 


gAo (in Coulomb gauge) 


Mv^ 


gA (in Coulomb gauge) 


Mv^ 



In (18), the balance between the kinetic energy and the potential energy rep- 
resented by the Aq term requires that gA^ scale like Mv^ . On the other hand, 
assuming that a gradient acting on Aq scales like Mw, (^^ requires that gAo 
scale like g'^Mv. These two estimates are consistent if the effective coupling 
constant as — .9^/47r at the scale Mv scales like v. This is identical to the naive 
estimate (|j) that followed from balancing the kinetic energy and the Coulomb 
term in the potential energy. Since this scaling relation follows simply from the 
consistency of the field equations, it applies to charmonium and bottomonium 
even though perturbation theory at the scale Mv is of questionable validity. 
The neglect of terms involving A in the field equations ([l^ ) and ([l9| ) is justified 
by the field equation for A, for which the dominant terms are 



(9o^ - V^) .gA 



gAoVgAo - ^V'^VV' ~ 0. 



(20) 



The last two terms in (po| ) scale like M^v^ and g'^Mv^, respectively, and they 
are comparable if g^ scales like v. Assuming that a gradient acting on A scales 
like Mv, we obtain the estimate Afw^ for gA in Table ^. 

Using the estimates for gA^ and gA in Coulomb gauge, we can obtain 
estimates for the field strengths. In Coulomb gauge, the dominant term in 
the chromoelectric field strength E is — VAq, and the resulting estimate for 
gE is M'^v^. The dominant term in the chromomagnetic field strength B is 
V X A, which leads to the estimate Af ^w** for gB. These estimates for gE and 
gH, although derived in Coulomb gauge, hold in general for matrix elements 
of gauge-invariant operators. 

According to the velocity-scaling rules in Table ^, the terms in the la- 
grangian density in (12) for minimal NRQCD scale like M*v^. Multiplying 
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by a volume factor of l/(Mw)'^, we find that quarkonium energies scale like 
Mv^. Each of the terms in 5C given in ( p^ ) scales like M^v^ and therefore 
contributes to onium energies at order Mv.KW other terms that can be added 
to the NRQCD lagrangian give contributions of order Mv^ or smaller. 

The validity of the velocity-scaling relations has been demonstrated con- 
vincingly by nonperturbative calculations of the bottomonium and charmo- 
nium spectrum using Monte Carlo simulations of lattice NRQCD Q The two 
parameters of minimal NRQCD can be tuned to give the spin-averaged spec- 
trum to an accuracy of about 30% for charmonium and about 10% for bot- 
tomonium. When the ^^-improvement terms are included, the errors decrease 
to about 10% for charmonium and to about 1% for bottomonium. These terms 
also give spin splittings that are accurate to about 30% for charmonium and 
to about 10% for bottomonium. 

2.5 Fock state expansion 

The simplest intuitive picture of quarkonium is that it is a bound state consist- 
ing of a Q and Q with very little probability of containing additional gluons or 
qq pairs. This simple picture is in fact realized in Coulomb gauge. Using the 



velocity-scaling rules of section 2.4, one can quantify the probabilities of Fock 
states containing additional gluons by determining how they scale with v. 

The Coulomb gauge (V • A = 0) is a physical gauge with no negative 
norm states, a necessary condition for a sensible Fock space. In this gauge, 
the scalar potential Aq does not propagate. Dynamical gluons are created and 



destroyed by the vector potential A. In Coulomb gauge, the lagrangian (12) 
can be reorganized as an expansion in powers of v. The powers of v can be 
made explici£l by rescaling the space-time coordinates r and t by l/{Mv) and 
l/(Afw^), respectively, rescaling the fields i(j and x by (Mu)'^/^, and rescaling 
the fields Aq and A by Mv^^^. The terms in the NRQCD lagrangian that are 
of order v^ after such a rescaling are 

Co = Ciight + V'^ (ido - gAo + ^ ) ^ + x^ Mo - gAo - ^ ) x- (21) 



2MJ " ^ y " ■:> ^' 2MJ 

This lagrangian in Coulomb gauge can be used to calculate quarkonium en- 
ergy levels to the same accuracy as the gauge-invariant lagrangian of minimal 
NRQCD. The terms in the NRQCD lagrangian that are of order v after the 
rescaling are 

+ charge conjugate terms. (22) 
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At first order in perturbation theory, these terms give transitions from 
Fock states to states that contain a dynamical giuon. The expectation values 
of the terms in ( |2^ ) between \QQ) Fock states vanish, so £i first contributes to 
quarkonium energy levels at second order in perturbation theory, giving shifts 
of order Mw'* . The terms in the Coulomb-gauge lagrangian that are of order 
v^ after rescaling are 



ijH-V'gAo)i^ - -J-i;^VgA„) x V • a^- 



+ ——'ip^igA X gA) ■ crip + charge conjugate terms. (23) 

The first term comes from expanding out the covariant derivative D in (12), 
while the last four terms in (p3) come from the ii^-improvement terms in (14). 
Energy levels calculated in Coulomb gauge using the lagrangian Cq + Ci + £2 
will differ only at order AIv^ from the energy levels calculated using the gauge- 
invariant NRQCD lagrangian (|2|) with 5C given by (pH). 

We now consider the Fock state expansion in Coulomb gauge for a quarko- 
nium state \H). It has the schematic form 

\H) -^^q\QQ) +^^qJQQ9) + ..., (24) 

where spin and color indices and momentum arguments have all been sup- 
pressed. The dominant Fock state \QQ) consists of a Q and Q in a color-singlet 
state with definite angular momentum quantum numbers '^ ~^^Lj. The higher 
Fock states, such as \QQg), include dynamical gluons or light qq pairs. Since 
the lagrangian Cq in ( pl| ) does not include any terms that couple A to V' or x, 
the probabilities for higher Fock states are suppressed by powers of v. 

The \QQg) states with the highest probabilities are those that couple to 
the dominant \QQ) state via the lagrangian £1 in (p^). We first consider 
the term ip''{igA ■ 'V)ip. We refer to a transition that proceeds via this term 
or its charge conjugate as an electric transition. An electric transition from 
the dominant \QQ) Fock state produces \QQg) states for which the angular- 
momentum quantum numbers of the QQ pair satisfy the selection rules AL — 
±1 and AS" = 0. The simplest way to determine the probabilities of these Fock 
states is to use the fact that a second-order perturbation in £1 changes the 
mass of the onium state by an amount of order Mw"*. This mass shift can be 
expressed as the product of the energy E of the virtual \QQg) state multiplied 
by its probability P. If the energy of the dynamical gluon is of order Mv, then 
E ~ Mv and we find that P ^ v^. If the gluon has energy of order Mv"^ or less, 
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then E ^ Mv^ and we obtain P ^ v^. We conclude that \QQg) states which 
satisfy the selection rules AL = ±1 and AS" — are dominated by very soft 
dynamical gluons with momenta of order Mw^ or less and have probabilities 
of order v'^. 

We next consider the term ■i/)'!' ( V x gA) -crilj in (|2^) . We refer to a transition 
that proceeds via this term or its charge conjugate as a magnetic transition. A 
magnetic transition from the dominant \QQ) Fock state produces \QQg) states 
that satisfy the selection rules AL = and AS = ±1. We can use the same 
argument as before to determine the probabilities of these Fock states, except 
that we must take into account the fact that the transition amplitude from the 
term tp^ (V x gA) • crtp is weighted by the momentum of the gluon. If the gluon 
has energy of order Mv, the mass shift from a second-order perturbation in 
£i is given correctly by the velocity-scaling rules to be of order Mv^. Since 
the virtual \QQg) state has energy E ^ Mv, we obtain a probability P ^ v^. 
The contribution to the mass shift from a gluon with energy of order Mv'^ is 
suppressed by a factor of w^ from the transition amplitudes and is therefore of 
order Mv^. Taking E ~ Mv^, we obtain P ^ v'^. We conclude that \QQg) 
states which satisfy the selection rules AL = and AS* = ±1 are dominated by 
dynamical gluons with momenta of order Mv and have probabilities of order 
v^. 

Similar arguments can be used to determine the magnitudes of the prob- 
abilities for other Fock states. Any such state can be reached by a sequence 
of electric transitions and zero or one magnetic transition. Electric transitions 
obey the selection rules AL = ±1 and AS* = 0, while magnetic transitions 
satisfy AL = and AS = ±1. Both electric and magnetic transitions change 
the color state of a color-singlet QQ pair to color-octet, and they change the 
color state of a color-octet QQ pair to either color-singlet or color-octet. The 
probability of a particular Fock state is determined by the color and angular- 
momentum quantum numbers of the QQ pair in that state. If that Fock state 
can be reached from the dominant \QQ) Fock state by a sequence of E electric 
transitions, then its probability scales like v'^^. If it can be reached by a se- 
quence of E electric transitions and a magnetic transition, then its probability 
scales like v'^^^^. 



2.6 Matching of NRQ CD and Q CD 

The NRQCD lagrangian contains adjustable parameters that must be tuned 
in order that its predictions for low-energy observables in the QQ sector agree 
with those of QCD. In the minimal NRQCD lagrangian, there are two param- 
eters: g and M. In addition to these parameters, the definition of NRQCD 
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requires an ultraviolet cutoff A to remove ultraviolet divergences. In the v^- 
improved lagrangian obtained by adding the terms in (|l4), there are 6 pa- 
rameters, (7, M, ci, C2, C3, and C4, in addition to the ultraviolet cutoff. The 
determination of the parameters in the NRQCD lagrangian is called matching. 

One could in principle determine the N parameters in the NRQCD la- 
grangian by tuning them so that the masses of N states in NRQCD match 
the corresponding masses in full QCD. For example, the parameters g and 
M of minimal NRQCD could be determined by matching the mass splittings 
between J/^ and "0' and between J/i/j and XcJ- Since the masses are sensitive 
to long-distance effects, they must be calculated nonperturbatively. The only 
reliable nonpcrturbativc method that is currently available is Monte Carlo sim- 
ulations of lattice NRQCD. While the masses in full QCD could in principle 
be computed nonperturbatively using lattice simulations, it is easier to take 
them directly from experiment. Using masses to tune the NRQCD parameters 
is an example of nonperturbative matching. 

Nonperturbative matching would become increasingly difficult as we strive 
for higher accuracy by adding more improvement terms. The determination 
of the parameters in the u^-improved lagrangian would require the nonper- 
turbative calculation of 6 masses as functions of 6 independent parameters. 
Fortunately, the asymptotic freedom of QCD provides an alternative, and that 
is perturbative matching. This matching procedure is based on the fact that 
QCD and NRQCD are equivalent except on distance scales of order 1/Af where 
perturbative QCD is by assumption accurate. The procedure for perturbative 
matching is the following: 

1. Use perturbative QCD to calculate scattering amplitudes between asymp- 
totic Q, Q, and gluon states with momenta k much less than M as func- 
tions of as and mg and expand them in powers of k/?TiQ. 

2. Use perturbative NRQCD to calculate the same scattering amplitudes in 
terms of the parameters in the NRQCD lagrangian and expand them in 
powers of k/Af . 

3. Adjust the NRQCD parameters so that the scattering amplitudes match 
to the desired order in k/mg, which we take to be of order v. 

It is essential to match scattering amplitudes or other physical observables 
rather than Green functions, because the construction of NRQCD involves 
field redefinitions. Such redefinitions can change the off-shell Green functions 
of the theory, but they leave on-shell physical observables unchanged. 

In present calculations in lattice NRQCD, the parameters are determined 
by a combination of nonperturbative and perturbative matching. The coeffi- 
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cients ci, C2, C3, and C4 of the w^-improvement terms are generally determined 
by perturbative matching, while the parameters g and M are determined by 
the nonperturbative matching of masses in the onium spectrum. Perturbative 
matching calculations can be used to relate these parameters to the fundamen- 
tal parameters of QCD. By combining these perturbative matching relations 
with lattice NRQCD calculations of the bottomonium spectrum, the QCD cou- 
pling constantj^s and the bottom quark mass rUb have been determined with 
high precision El 

The method of perturbative matching is somewhat paradoxical. We have 
assumed that M is large enough that perturbation theory is accurate at the 
scale M . We allow for the scale Mv to be small enough that perturbation 
theory is not reliable at that scale. If that is the case, perturbative calcula- 
tions in NRQCD would never give accurate results for physical observables, 
since NRQCD only reproduces full QCD accurately at scales of order Mv or 
less. Nevertheless, a comparison of perturbative calculations in NRQCD and 
full QCD can be used to accurately determine the parameters in the NRQCD 
lagrangian. The reason for this is that the tuning of the parameters of NRQCD 
that makes this theory equivalent to QCD at momenta of order Mv or smaller 
also makes the perturbative approximations to these theories equivalent. Per- 
turbation theory breaks down in precisely the same way for both theories, 
predicting among other things, the existence of asymptotic states consisting of 
isolated quarks and gluons. Since the parameters in the NRQCD lagrangian 
are sensitive only to momenta on the order of M where perturbative QCD is 
accurate, they can be correctly determined by matching perturbative calcula- 
tions in QCD and NRQCD. 

As an illustration of perturbative matching, we consider the simplest per- 
turbative observable. This is the energy-momentum relation for the heavy 
quark, which is given by the location of the pole in the heavy-quark propaga- 
tor. At tree level in full QCD, the energy-momentum relation is 



E^f^^^^.n^^^£^^^.... (25) 



At tree level in NRQCD, we can read off the energy momentum relation from 
the lagrangian ([12|): 

2 4 

E = ^ ci^-- + .... (26) 

2Af ^8M3 ^ ' 



By matching the expressions ( |25| ) and (26) we find 

M = mq, ci = 1. (27) 
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If the energy-momentum relations are computed to higher order in perturba- 
tion theory, the matching will give perturbative corrections to the results in 
(P7|). Since the parameters M and ci are sensitive only to short distances of 
order l/mg or smaller, the corrections can be expressed as power series in 

3 Annihilation Decays of Heavy Quarkonium 

3. 1 Decay of rjc in the Color-Singlet Model 

A simple intuitive picture of the annihilation decay of a quarkonium state is 
that it proceeds through the annihilation of the QQ pair in the dominant Fock 
state into gluons and light qq pairs. These light partons ultimately hadronize 
into the observed final states that consist of light hadrons. The inclusive an- 
nihilation rate of the QQ pair can be plausibly calculated using perturbative 
QCD. By combining that perturbative calculation with a phenomenological 
wavefunction for the dominant \QQ) Fock state, we can calculate the anni- 
hilation decay rate of the quarkonium. This model for annihilation decays is 
called the color-singlet model. 

The simplest illustration of the color-singlet model is the calculation of 
the decay rate of the rjc- In the color-singlet model, the rjc is modeled by a 
cc pair in a color-singlet ^Sq state. Its wavefunction is the product of a color 
factor 5ij/^/?>, a spin factor (ti — it)/'\/2; and a coordinate-space wavefunction 
■^(r) — R(r)/^f^. A cc pair in such a state can annihilate into two gluons. If 
we assume that the two gluons hadronize into light hadrons with probability 
1, the decay rate can be written as 

The T-matrix element for this decay can be expressed in terms of the momentum- 
space wavefunction V'(q) of the rjc- 



1 /■ d^q 

^2mZj W? 



T[,7, ^ 5(k)g(-k)] = ^== / ^ ^(q) T[c(q)c-(-q) -> .g(k)g(-k)], 



(29) 

where we have suppressed all color and spin indices. The cc annihilation am- 
plitude T varies significantly with q only when |q| is on the order of rric or 
larger. The wavefunction has significant support only for |q| of order rricV. 
Since T is almost independent of q for such small values of |q|, we can set 
q = in the annihilation amplitude. The resulting expression for the decay 
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rate has a factored form: 






.3, 2.SiM,^-2l^^^^^^^_^^^^^^^^ .^^^^^ 



(27r)32fc Mr,, 

The integral over q gives the wavefunction evaluated at the origin, i?(0)/\/47r- 
The integral over k in (^0|) can be calculated from the lowest order QCD 
Feynman diagrams for cc —f gg. The final expression for the decay rate is 

r(,j » ^mo)f^ (31) 

3.2 Decay of t]c m NRQCD 

The decay rate of the rjc into individual final states consisting of light hadrons 
can not be described within the framework of NRQCD. One obstacle is that 
we have imposed a symmetry on the effective field theory that guarantees the 
separate conservation of the numbers of c's and c's and therefore forbids the 
annihilation process cc — > gg. Furthermore, in the construction of NRQCD, 
we have integrated out gluons with momenta on the order of rric. Even if we 
relax the definition of the effective theory to allow gluons with momenta of 
order rric and interaction terms in the lagrangian that allow cc annihilation, 
we cannot described annihilation decays accurately, because the interaction of 
a c or c with a gluon of momentum rric cannot be described accurately in a 
local nonrelativistic theory. 

While the decay rate of the r]c into a specific final state consisting of light 
hadrons cannot be described within NRQCD, the total decay rate can. The 
conservation of the number of c and c is not an obstacle, because the optical 
theorem can be used to express the inclusive annihilation rate in terms of an 
amplitude that conserves the numbers of c and c. As an illustration, applying 
the optical theorem to the expression (|2q), we find that the decay rate can be 
written 



r(^c 



f f (Pq f <fq' 



(2M^J2 J (27r)3 J (27r)3 

xV'(q) 2Imr[c(q)c(-q) ^ c(q')c(-q')] ^*(q')- (32) 



The process cc —f cc conserves the numbers of c's and c's. The factored form 
( pO| ) is recovered by using the fact that the imaginary part of the T-matrix 
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element for cc —^ cc is insensitive to momenta q and q' on the order of Mv 
where the wavefunctions have their support: 



nvc) ^ 



(2M^e)' 






2 ImT [c(0)c(0) -^ c(0)c(0)] . (33) 



The integral over q in (|3^) gives the square of the wavefunction evalu- 
ated at the origin. This reflects the fact that the initial c and c must have 
spacetime separations of order 1/M in order to annihilate, and this separation 
is small compared to the length scale \/{Mv) of the wavefunction. That the 
cc pair must have a spacetime separation of order \/M follows from the fact 
that the Feynman diagrams for the annihilation process involve a heavy quark 
propagator that is off its mass-shell by an amount of order M . The same ar- 
gument implies that the final c and c in the T-matrix element in (BS) must 
have a space-time separation of order \/M . While it is not quite as obvious, 
the space-time separation of the annihilation points for the initial cc pair and 
the final cc pair must also be of order 1/M. This follows from the requirement 
that the wavefunctions of the annihilation gluons must overlap, which localizes 
the production point of a gluon to within its wavelength 1/M. 

The fact that cc annihilation occurs within a region whose size is of order 
1/M provides a clue as to how the effects of annihilation can be taken into 
account in NRQCD. All modes with momenta of order M that can be sensitive 
to the length scale 1/M have been removed from this effective theory. Thus 
the effects of the annihilation can be reproduced by including in the NRQCD 
lagrangian (n2|) a local 4-fermion interaction term that destroys a cc pair and 
creates it again. The specific term that is relevant to rjc decay is 

5C = -^i^^xxH. (34) 



where the coefficient / is dimensionless. The term (34) annihilates a cc pair in 
a color-singlet ^5*0 state and then creates a cc pair in the same state. 

The dimensionless coefficient / in ( p4[) can be determined by perturbative 
matching of the cc -^ cc scattering amplitudes in full QCD and NRQCD. 
In full QCD, the scattering amplitude includes box diagrams of order a1 in 
which the scattering proceeds through intermediate states consisting of two 
gluons. In NRQCD, this contribution to the scattering amplitude can only be 
reproduced by 4-fermion interactions such as those in (p4|). By matching the 
NRQCD scattering amplitude from the term (|3j) with the annihilation part 
of the scattering amplitude in full QCD, we can determine the coefficient /. 
Since the c and c can annihilate into two on-shell gluons, the QCD scattering 
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amplitude has an imaginary part. The coefBcient / in (p4) must therefore have 
an imaginary part and it is particularly simple to calculate. The result is 

Im/ = ^^^. (35) 

Since this coefficient is sensitive only to distances of order 1/mc or smaller, 
the running coupling constant is evaluated at the scale rric. The fact that 
coefRcients in the NRQCD lagrangian have imaginary parts implies that the 
hamiltonian for NRQCD is not hcrmitian. This is perfectly natural, since we 
have removed states from the theory that are essential for exact unitarity. In 
particular, we have eliminated the light partons with momenta on the order of 
rric that can be produced by the annihilation of the c and c. 

We now consider the effect of the correction terms ( |34| ) on the energy of the 
•qc- If that term is treated as a first-order perturbation, the resulting correction 
to the energy of the rjc is 

^^"= ^ -jp — 2m;^ — ' ^^^) 

where we have assumed that the state \rjc) has the standard relativistic nor- 
malization. Since the coefficient / has an imaginary part, this energy shift has 
an imaginary part. A state whose energy has a small imaginary part —T/2 
should be interpreted as a resonance of width F. Thus the width of the rjc due 
to the term ( p^ ) is 

r(,e) = ^]^^%^(^c|^^xxVhc). (37) 



The connection with the result (|3^) from the color-singlet model is made by 
inserting a complete set of states between x ^tnd x^ in the matrix element. 
Assuming that the sum over states is dominated by the vacuum, we have 

{ric\^^XxH\ric) ~ \{()\xH\Vc)\^ ■ (38) 

In the color-singlet model, the yyc-to-vacuum matrix element on the right side 
of (^) can be expressed in terms of the wavefunction: 

3" 



(OlxVhc) « ./2MZ\J^R{0). (39) 

Inserting (^) into (^^ and identifying M^^ with 2mc, we reproduce the result 
( p^ ) from the color-singlet model. 
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3.3 NRQCD Factorization Formula 

In the expression ( p7| ) for the decay rate of the rye, short-distance and long- 
distance effects have been factored. Long-distance effects involving the quarko- 
nium wavefunction appear only in the NRQCD matrix element, which scales 
like M^v'^ according to the velocity-scaling rules in Table 0. Short-distance 
effects involving the annihilation of the QQ pair appear only in the coefficient 
2Im//M^, which is expressed in terms of the fundamental parameters as and 
771c of QCD. Thus the expression for the decay rate in (pTJ) scales like alv'^ . 

The formula ( p7| ) can be generalized to all orders in a^ and to all orders 
in V. The general factorization formula for the annihilation decay rate of a 
quarkonium state H is 

nH) - l^Y.^ran{H\Omn\H), (40) 

ran 

where the onium state \H) — \H{V — 0)) has the standard relativistic nor- 
malization. The sum in (BO) extends over all operators that can appear in the 
NRQCD lagrangian and that have the form 

Omn = V'^^rnXX^^nV'- (41) 

These operators must be gauge invariant, invariant under parity and charge- 
conjugation, and scalars under rotations. Each of the factors /C„ and /Cm is the 
product of a spin-matrix (1 or cr*), a color matrix (1 or T"), and a polynomial 
in D and [Dq, D] — 5E. The operator Omn in ( pl| ) annihilates a QQ pair in a 
color and angular-momentum state determined by /C„ and creates a QQ pair 
at the same point in a state determined by /Cm. 

The NRQCD factorization formula ( pO| ) untangles the effects of short dis- 
tances of order \/M from those of long distances of order \/{Mv) or larger. 
All long-distance effects involving the quarkonium wavefunction are factored 
into the NRQCD matrix elements. Short-distance effects involving the annihi- 
lation of the QQ pair are contained in the coefficients. The coefficient Cmn in 
( pO| ) is twice the imaginary part of the coefficient of the operator Omn in the 
NRQCD lagrangian. If that operator has scaling dimension dmm then Cmn is 
l/mj""" multiplied by a power series in as{'mQ). 

The coefficients Cmn can be calculated using perturbative matching meth- 
ods. A general matching prescription, called the threshold expansion method, 
has been developed by Braaten and ChenQ The matching calculations are car- 
ried out using perturbative asymptotic states cc = cc(q, ^, rj) that consist of a c 
and a c with relative momentum q and in a spin/color state that is represented 
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by the Pauli spinors ^ and rj. The matching prescription is 



X 



pQCD 



■"^-^ pNRQCD 



where P = {2^Jm'^ + q^, 0) is the four-momentum of the cc pair and T^c^x 
is the T-matrix element for its annihilation into a final state X consisting of 
light partons. The sum over X on the left side of ( [42| ) includes integration 
over the phase space of the light partons. The complete determination of the 
short-distance coefficients requires the use of different states cc and cc! in the 
T-matrix element and in its complex conjugate. In the matching procedure, 



the left side of (42) is calculated using perturbative QCD, and then expanded 
in powers of the relative momenta q and q'. The matrix elements on the right 
side are calculated using perturbative NRQCD, and then expanded in powers 
of q and q'. The coefficients Cmn are determined by matching these expansions 
order by order in as- 

The relative importance of the various terms in the factorization formula 
( p^ ) is determined by the order in as of the coefficient Cmn and by the order 
in V of the matrix element. The magnitude of the matrix element can be 
estimated using the velocity-scaling rules for operators given in Table and 
the estimates for the probabilities of higher Fock states in Coulomb gauge 



given in Section 2.5. If the operator Omn annihilates and creates a QQ pair 
in the same color and angular- momentum state as in the dominant \QQ) Fock 
state of iJ, then the estimate for the matrix element is obtained by dividing the 
estimate for the operator Omn from Table y by M^w^. Otherwise, we must take 
into account suppression factors from the transitions required to go from the 
dominant \QQ) state to a Fock state in which the QQ pair can be annihilated 
by the operator Omn and then back to the dominant \QQ) state. There is a 
suppression factor of v for every electric transition and a suppression factor 
of v'^/'^ for every magnetic transition. Let E and M be the total number of 
electric and magnetic transitions required. If /C„i and /C„ contain D factors of 
the covariant derivative D but no factors of (^E or gB, then the matrix element 
{H\Omn\H) scales like «3+D+B+3Af/2^ Yov each factor of ^E or 5B, there is an 
additional suppression factor of v'^ or w^, respectively. 

Spin symmetry relates NRQCD matrix elements for quarkonium states 
that differ only in their spin quantum numbers. An example involving the rjc 
and J /ip (the lowest ^5*0 and '^ Si states of charmonium) is 

{J/i^\i>^ax-X^cTi^\J/^) ~ {Tlc\i>^XxH\Tlc). (43) 
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Since spin symmetry is only an approximate symmetry of NRQCD that is 
broken at order v^, the equality (E^ holds only up to corrections of relative 
order u^. 

The vacuum-saturation approximation can be used to express some of 
the matrix elements in the factorization formula ( [40| ) in terms of vacuum-to- 
quarkonium matrix elements. This approximation can only be applied if the 
operator Omn annihilates and creates QQ pairs in the same color and angular- 
momentum state as in the dominant \QQ) Fock state of H. In this case, we 
can insert a complete set of states between x f^nd x' ■ 

{H\ij^lC,nXX^IC^^\H) = ^(i/lV-^/C^xl^X^lx^-^nV^I^)- (44) 

X 

The vacuum-saturation approximation consists of keeping only the vacuum 
term |0) (0| in the sum over states, as illustrated in (|38|). The vacuum-saturation 
approximation is a controlled approximation with an error of relative order i;^. 
This follows from the fact that in Coulomb gauge, the next most important 
term in the sum over states in (|44| ) is a \gg) Fock state, which contains two 



dynamical gluons. The leading contributions to the matrix elements in (44) 
then come from \QQgg) Fock states whose probabilities are of order w^. 

The NRQCD matrix elements are sensitive to long-distance effects, and 
therefore can only be calculated using nonperturbative methods. The only 
practical nonperturbative method that is presently available is Monte Carlo 
simulations of lattice NRQCD. The first such calculations have been car- 
ried out recently by Bodwin, Sinclair, and KimB They demonstrated that 
the relation ( pq ) implied by the vacuum-saturation approximation holds to 
within numerical accuracy. They calculated the matrix elements (0|x^V'l'7c) 
and {Q\x^T)'^'4>\rjc) and their analogues for bottomonium. These matrix ele- 
ments contribute to decays of the lowest S-wave states at leading order in v and 
at relative order v^ . They also calculated (0|x^D • cnp\xcQ) and (xco|'0^'''2^°X ' 
X^crT^f/^IXco) and the analogues of these matrix elements for bottomonium. 



As discussed in Section 3.4, these are the matrix elements that contribute to 
decays of the lowest P-wave states at leading order in v. Thus far, the cal- 
culations of NRQCD matrix elements have been carried out only in minimal 
NRQCD and without dynamical quarks. 

3.4 Annihilation Decays 0/ P-wave States 

The NRQCD factorization formula for annihilation decay rates has dramatic 
implications for the decays of P-wave states, such as he (the ^Pi state of 
charmonium) and XcJ, J = 0, 1, 2 (the ^Pj states). Calculations of their anni- 
hilation decay rates in the color-singlet model suffer from infrared divergences. 

23 



The NRQCD factorization approach not only resolves the problem of the in- 
frared divergences, but it also leads to new qualitative insights about P-wave 
charnionium. 

We first consider the Fock state expansion of he and Xcj in Coulomb gauge. 
The dominant |cc) Fock state consists of a color-singlet cc pair in a -'^ Pi state 
for he and a '^Pj state for XcJ- An electric transition from the dominant Fock 
state produces a \ccg) state, with the cc pair in a color-octet state with angular 
momentum quantum numbers ^5*0 or ^D2 for he and ^Si, '^-Di, ^D2, or ^D^ for 
XcJ- Therefore these Fock states have probabilities of order w^. A magnetic 
transition from the dominant \cc) Fock state produces a \ccg) state, with the 
cc pair in a color-octet state with angular-momentum quantum numbers '^Pj 
for he and ^Pi for Xcj- This Fock state has a probability of order v'^. All other 
Fock states have probabilities of order v'^ or smaller. 

We proceed to identify the most important matrix elements in the NRQCD 
factorization formula. The lowest dimension operator tjiat can annihilate the 
cc pair in the dominant Fock state of the he is ■(/'^DX'X^D-0, where D is defined 
by x^D'i/) = xH'^i^) ~ (I^x)^V'- For XcJ, the lowest dimension operator has the 
form ■0'''D'"cr"xx'''DV'''0, where the indices are contracted in different ways for 
J — 0,1, and 2. The matrix elements of these operators for /icXcOiXcij E^nd 
Xc2 are related by spin symmetry. Up to corrections of relative order v'^ , they 
satisfy 

{he\ijH-^b)x-xH-^i>)i^\hc) 

~ i(Xco|^^(-§D. ^)xxH-|D. ^)^|Xco) 

~ ^{XciHH-i^ X ct)x ■ xt(-§D X <t)7Mxci) 

« (Xc2|^^(-|5)("a")xxt(-i5)("V")VIXc2), (45) 

where pC™") denotes the symmetric traceless part of a tensor p™". The 
vacuum-saturation approximation can be used to express these matrix ele- 
ment in a simpler form. Up to corrections of relative order w^, the matrix 
element for the Xco can be written 

(Xco|V^^(-tD- <t)xx^(-|D. a)^AXco) « |(0|xt(-^D • <T)V|Xco)f • (46) 

In the color-singlet model, the vacuum-to-XcO matrix element on the right 
side of (^) can be expressed in terms of the radial wavefunction R{r) for the 
P-wave states: 

(Olxn-tD.^)^IXco) « v/2M;^y^P'(0). (47) 
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According to the velocity-scaling rules, the matrix elements in (E^) scale 
like v^ . For a consistent analysis, we must also include all other matrix elements 
that scale like v^. By enumerating the possibilities, one can see that the only 
other operators whose matrix elements scale like v^ are ip'^T'^xX^T'^i' for he 
and ip^ (tT°'x ' X '^T°"4) for XcJ- These operators annihilate and create cc pairs 
in color-octet ^5*0 and ^Si states, respectively. Spin symmetry implies that the 
matrix elements for h^, Xco, Xci, and X2 are equal, up to corrections of relative 
order v^: 

{K\i^^T-XX^T'^^\hc) ~ {xcA^^^crT'^X-X^tTT'^i'Wcj), J = 0,1, 2. (48) 

We have found that, up to corrections that are suppressed by u^, the 
annihilation decay rates of the P-wave states can all be expressed in terms of 
the following two independent matrix elements: 

(0i> ^ KT-r > (49) 



2Mxco 

2-'"XcO 

Their short-distance coefficients can be calculated as power series in as{'mc)- 
The annihilation processes that contribute to the coefficients at order a1 are 
cc — > gg and cc — > qq. If we keep only those terms in the short-distance 
coefficients, the annihilation decay rate are 

r(xco) - '^^^m + ^-^^^^m , (52) 



3m^ 

r(xci) - "'TT' m. (53) 



nfiTa1{mc) 
3'm'?. 



r(Xc.) « ^^^^(C)r) + ^^^^%^(08), (54) 

where ny = 3 is the number of flavors of light quarks. At this order in as, 
the decay rates for he and Xci receive contributions from (Os) only, because 
Yang's theorem forbids the annihilation process cc — > gg for a cc pair in a state 
with total angular momentum 1. At order a^, all of the P-wave states have 
contributions from (Oi). At this order in a^, the short-distance coefficients of 
(Oi) depend logarithmically on a factorization scale /z that can be interpreted 
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as an infrared cutoff on the energy of soft gluons. The matrix element (Og) 
also depends logarithmically on ^, which in this case can be identified with the 
ultraviolet cutoff ofJNRQCD. The /x-dependence cancels between (Os) and the 
coefficient of (Oi)E£lln the color-singlet model, the (Og) terms are absent and 
the decay rate depends logarithmically on the infrared cutoff /i. The NRQCD 
factorization approach provides a simple and natural solution to this problem. 
The NRQCD factorization formula for a P-wave state has a simple in- 
terpretation in terms of the Fock state expansion in Coulomb gauge. The 
color-singlet terms proportional to (Oi) are contributions from the dominant 
|cc) Fock state, while the color-octet terms proportional to {0%) are contribu- 
tions from a \ccg) Fock state. The \ccg) state has a small probability of order 
v^ , and its effects on most observables are small compared to those of the | cc) 
state. However, in the case of the annihilation decay rate, the effects of the |cc) 
state are suppressed by u^ due to the orbital angular momentum of the |cc) 
pair. The contribution of the \ccg) Fock state has no angular-momentum sup- 
pression and therefore contributes at the same order in v. For the he and Xco, 
the effects of the \ccg) Fock state actually dominate, because the short-distance 
coefficients of (Oi) are suppressed by a factor of as{mc). 

4 Inclusive Production of Heavy Quarkonium 

4-1 Topological Factorization 

The cross section for producing a quarkonium state iJ in a high energy process 
necessarily involves both "short distances" of order 1/M or smaller and "long 
distances" of order l/(Mv) or larger. The creation of the QQ pair involves 
short distances, because the parton processes that produce the QQ pair always 
involve particles that are off their mass shells by amounts of order M and can 
therefore propagate only over short distances. The binding of the Q and Q 
into the state H involves long distances, because gluons whose wavelengths 
are comparable to or larger than the size of the bound state, which is of order 
l/{Mv), play a large role in the binding. 

The production of quarkonium in a high energy physics process typically 
involves another hard momentum scale Q in addition to the scale M. If the 
production cross section is sufficiently inclusive, it can be described by an 
NRQCD factorization formula that separates short-distance effects involving 
the momentum scales Q and AI from long-distance effects that involve lower 
momentum scales. It is convenient to separate the derivation of the factoriza- 
tion formula into two steps. In the first step, which we refer to as topological 
factorization, the standard factorization methods of perturbative QCD are 
used to separate the effects of the hard momentum scale Q from those of the 

26 



soft momentum scale Aqcd- These methods are applicable even if we identify 
the hard scale Q with the heavy quark mass M . An additional step is required 
to separate the effects of the scale M from those of the scale Mv. 

The general expression for the inclusive cross section for the production of 
a quarkonium state H with four-momentum P is 

V da(12 -^ H(P) +X) = ^- , ^^f 

^ ^ "■ > > AE1E2V12 {2TTf2Ep 

X Y, (27r)4j4(fci +k2-P-kx) \T,2^HiP)+x\^ , (55) 

X 

where Ti2^h{p)+x is a T-matrix element for producing H and the additional 
particles X and the sum on the right side includes integration over the phase 
space of the additional particles. At the parton level, both Ti2^h(p)+x and its 
complex conjugate can be expressed as sums of Feynman diagrams. The prod- 
uct of a single diagram in 7^2^ff(p)-i-x and a single diagram in 'T*^^^,p.,-^ 
is called a "cut Feynman diagram." Using the factorization methods of per- 
turbative QCD, one can identify the cut diagrams that dominate in the limit 
Q ^ 00. After taking into account cancellations between real and virtual soft 
gluons, the dominant cut diagrams have the following structure: 

• a hard- scattering subdiagram Ti to the left of the cut. The outgoing 
lines include a QQ pair with small relative momentum of order Mv and 
additional hard partons. There can also be incoming hard-parton lines 
if the process involves hadrons in the initial state. 

• a hard- scattering subdiagram Ti* to the right of the cut that is just the 
mirror image of 7i. 

• a jet-like subdiagram Ji for each of the hard partons attached to W. The 
subdiagram extends through the cut and is attached to Ti. and to Ti* by 
single hard parton lines. 

• an oniuni subdiagram O that extends through the cut and is attached to 
H and to Ti* by Q and Q lines that have small relative momentum. 

The cut diagrams that do not have the above structure are suppressed by pow- 
ers of 1/Q. With this topological factorization of the dominant cut diagrams, 
all effects involving the hard momentum scale Q are factored into the hard- 
scattering subdiagrams H and H* , while all effects of the soft scale Aqcd are 
factored into O, J7i, J2, ■ ■ ■■ The gluon interactions that bind the QQ pair 
into the onium state H are also contained within the onium subdiagram O. 
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The proofs of the factorization theorems of perturbative QCD are very 
difhcult, and exphcit proofs are available only for a very few processes, such 
as inclusive hadron production in e~^e~ annihilaiion and the Drell-Yan process 
for lepton pair production in hadron coUisionsLJ However there is no apparent 
obstacle to extending these proofs to inclusive onium production. 

4.2 NRQCD Factorization 

After topological factorization, the effects of the scale Mv are distributed in a 
complicated way between 7i, 7i*, O, and the Q and Q propagators that connect 
them. The onium subdiagram O involves the scale Mv because gluons with 
momentum of order Mv play an important role in the binding of the QQ pair 
into an onium state. The hard-scattering subdiagrams Ti and Ti* involve the 
scale Mv, because the outgoing Q and Q lines have relative momenta on the 
order of Mv. It is in factoring the scale Mv out of Ti and Ti* that NRQCD 
enters into the picture. 

Consider the part of the cut diagram that includes 7i, O, and the Q 
and Q propagators that connect them. If the Q and Q have four-momenta 
^P + q and ^P~Q7 then the diagram also involves an integral over the relative 
momentum q. A simple way to disentangle the momentum scale Mv from Ti 
is to expand it as a Taylor series in q and absorb the factors of q as well as 
the integration over q into O. Each term in the Taylor expansion corresponds 
to a local operator that creates a QQ pair from the vacuum. By applying a 
similar procedure to disentangle the momentum scale Mv from Ti* , the onium 
subdiagram is reduced to vacuum-expectation values of local operators that 
create and annihilate QQ pairs. After a renormalization group transformation 
and appropriate field redefinitions, the matrix elements can be expressed as 
expectation values in the NRQCD vacuum of the form 

{O^n) = {0\x^JC„,^pVH^^IC.nx\0), (56) 

where Vh projects onto states that in the asymptotic future contain the 
quarkonium state H plus soft partons S whose total energy is less than the 
ultraviolet cutoff of NRQCD: 

Vh - ^|ff + 5)(if + 5|. (57) 



If we integrate the dominant cut diagrams over the phase space of all the 
hard partons in the process, we obtain the NRQCD factorization formula for 
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Ml2^H{P)+X) = TTT^jr— 7:t-p^E^-(^1'^^'^)(^™. 



the inclusive cross section: 

(58) 
The sum in (BSl) extends over aU NRQCD matrix elements of the form (pq). 
The product of the operators x^^mV' and ^^ICnX must be gauge-invariant. 
It need not be rotationally invariant if the quarkonium state H is polarized. 
In the factorization formula (pq), the hard-scattering subdiagrams Ti and Ti* 
and the jet-like subdiagrams J7i have all been subsumed in the short-distance 
coefficients Cmn- 

Since the coefficients Cmn in ( |58| ) involve only short distances of order 
l/M or larger, they can be expressed as perturbation series in as{M). The 
threshold expansion methcid provides a general prescription for calculating the 
short-distance coefficients^ Denoting by cc{P) a state consisting of a c and c 
with relative momentum q that has been boosted to four-momentum P, the 
matching prescription is 



E {2^f5\k, +k2-P~-kx) {T,^^c-c-iP)+xTr, 



cc'(P)+X) ^l2^cc(P)+X 



X 



pQCD 



~ yCmnikl,k2,P)(0\x^lCmi^Vcc',cc^^lCnm , (59) 

^^ pNRQCD 

ran 

where the projection operator in the NRQCD matrix element is 

Pcc'.ce = E |cc' + 5) {cc + S\ . (60) 



The left side of (^9|) is to be calculated using perturbative QCD, and then ex- 
panded in powers of the relative momenta q and q' of the cc pairs. The matrix 
elements on the right side are to be calculated using perturbative NRQCD, 
and then expanded in powers of q and q'. The coefficients Cmn are then 
determined by matching these expansions order by order in as. 

The relative importance of the various terms in the factorization formula 
( [lO|) is determined by the magnitudes of the coefficients Cmn and by the order in 
V of the matrix elements (O^^). The size of the coefficient Cmn is determined 
not only by the order in as , but also by its dependence on dimensionless ratios 
of kinematic variables that are involved in the cc production process. The 
magnitudes of the matrix elements can be estimated by using the velocity- 
scaling rules for operators given in Table and the scaling with v of the rates for 
electric and magnetic transitions. If the operator Omn creates and annihilates 
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a QQ pair in the same color and angular-niomcntum state as in the dominant 
\QQ) Fock state of H, then the magnitude of the matrix element (Om„) is 
estimated by multiplying the factors in Table || and dividing by M'^v^. For 
other matrix elements, we must take into account suppression factors from the 
transitions required to go from the QQ state created by the operator -0t/C„x 
to a state in which the QQ pair has the same quantum numbers as in the 
dominant \QQ) Fock state and then to a state in which the QQ pair can be 
annihilated by the operator x t^mip- There is a suppression factor of v for every 
electric transition that is required and a suppression factor of v^'"^ for every 
magnetic transition. The scaling of the production matrix element (Omn) with 
V is identical to that of the corresponding decay matrix element {H\Omn\H). 
The NRQCD matrix elements that appear in the NRQCD factorization 
formula (g8|) can be simplified by using symmetries of NRQCD. Rotational 
symmetry is an exact symmetry of NRQCD. It implies, for example, that 

ix^a^T^^Vj/^^^a'T'^x) = ^(x^^T^ T',//^ ^V^T'^x)- (61) 

Spin symmetry is an approximate symmetry of NRQCD that holds up to cor- 
rections of order v'^. It implies, for example, that 

ix^a^T'^ijr^^^^ij^a'T-x) ~ UxjU^ix^T'^ib V,,^ ^^T'^x), (62) 

where A specifies the polarization of the J/ip and Uix is the unitary matrix 
that transforms vectors from the spherical basis to the Cartesian basis. 

The vacuum-saturation approximation can be used to simplify the matrix 
elements of operators that create and annihilate cc pairs in the dominant Fock 
state of the quarkonium. In the vacuum-saturation approximation, the projec- 
tion operator Ph defined in ( p7| ) is replaced by the single term \H)(H\. This 
is a controlled approximation in NRQCD, holding up to corrections that are 
of order v^. It implies, for example, that 

(OlxVV ^^(A) ^^'T^XIO) « (0|xVV|^(A))(V^(A)|VV^x|0). (63) 

The matrix elements in the NRQCD factorization formula involve long- 
distance effects and therefore can only be calculated using nonperturbative 
methods. Unfortunately, in contrast to the decay matrix elements {H\0„in\H), 
there are no effective prescriptions for calculating the production matrix ele- 
ments (O^n) using lattice NRQCD. The problem lies in implementing on the 
lattice the projection defined by {pl\). Thus these NRQCD matrix elements 
must be treated as phenomenological parameters to be determined by exper- 
iment. The only exceptions are the matrix elements to which the vacuum- 
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saturation approximation can be applied. Vacuum-to-quarkonium matrix ele- 
ments of the form {Hlip'' ICnX\0) can be calculated using Monte-Carlo simula- 
tions of NRQCD. 

4.3 Prompt Charmonium at the Tevatron 



The NRQCD factorization framework i58) has many applications, some of 
which are described in a recent reviewEHOne application for which the im- 
plications are particularly dramatic is the production of prompt charmonium 
at large transverse momentum in pp collisions. At sufficiently large trans- 
verse momentum pT, the cross section for pp —f ip -\- X is dominated by gluon 
fragmentationE3 It can be factored into the cross section for producing a gluon 
with transverse momentum pr/z and a fragmentation function Dg^^{z) that 
gives the probability that the jet initiated by the gluon includes a tA carrying 
a fraction z of the gluon momentum: 

da{pp ^ iP{P) + X) = f dzda{pP^g{P/z)+X)Dg^^{z). (64) 

Jo 

Using the NRQCD factorization approach, the fragmentation function can be 
expressed in the form 

Dg^^{z) - 5]d„„(z)(OL), (65) 

van 

where all effects of the momentum scale rric have been factored into the short- 
distance coefficients dmn{z)- The relative importance of the various terms 
in the fragmentation function is determined by the order in v of the matrix 
element and the order in as of its coefficient. 

The matrix element that is leading order in v is {x^a^ijjV^ij'^a^x), which 
scales like v^ . This term corresponds to the formation of a t/j from a cc pair that 
is created in a color-singlet ^5'i state, and it is the only contribution in the color- 
singlet model. The leading contribution to the short-distance coefficient of this 
term in the gluon fragmentation function ( p5|) is of order a^ and comes from 
the parton process g* — > ccgg. Keeping only this term in the fragmentation 
function (|5|), the cross section predicted by ( |64|) is about a factor of 30 below 
recent data on prompt ip production at the Tevatron. 

The color-singlet-model term in the gluon fragmentation function scales 
like a^f'^. All other terms have matrix elements that scale like v^ or smaller. 
There are however terms whose short-distance coefficients are suppressed by 
fewer powers of as ■ There is one coefficient in particular that is of order Us , 
because it receives a contribution from the parton process g* -^ cc. The matrix 
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element is {x' a''T°'ip'P^ip'' a''T"-x) and it scales like w^. This term corresponds 
to the formation of a V' from a cc pair that is created in a color-octet "^^i state. 
At leading order in a^, this term in the fragmentation function is 

D,_,4z) = ^^^^S{1 - z) (xV^T"^ V^ i'^a'^T^x) ■ (66) 

Braaten and Fleming proposed that the enhancement from the two fewer pow- 
ers of as relative to the color-singlet model term can overcome the suppressian 
by w'^, and that this term might dominate the gluon fragmentation functionEj 
The pT-dependence predicted by this mechanism is in agreement with the 
Tevatron data. The normalization depends on the unknown matrix element in 
( pq ) , but the value of the matrix element required to fit the data is consistent 
with suppression by a factor of v'^ relative to the color-singlet-model matrix 
element {x^a^ij^V^iji'^cr^x)- 

Cho and Wise pointed out that this ppSjduction mechanism has dramatic 
implications for the polarization of the ■0" At leading order in a^, the ^'s 
produced by gluon fragmentation will be 100% transversely polarized. The 
radiative corrections to the fragmentration function were examined by Beneke 
and Rothstein, and they jConcluded that the spin alignment at large pt will 
remain greater than 90%ll3 The largest corrections to the spin alignment at 
values of Pt that can be measured at the Tevatron come from nonfragmentation 
contributions that fall like l/p^j P^d these contributions have recently been 
calculated by Beneke and KraemertZI An experimental measurement of the spin 
alignment in agreement with these predictions would constitute a dramatic 
triumph of the NRQCD factorization approach. 

5 Conclusions 



The NRQCD factorization formulas (40) and ( p8| ) provide a firm theoreti 



cal foundation for analyzing annihilation decay rates and inclusive production 
rates of heavy quarkonium. The short-distance coefficients can be calculated 
as power series in the running coupling constant as(Af ), and the long-distance 
factors are defined in terms of NRQCD matrix elements that scale in a definite 
way with v. This approach not only provides a framework for carrying out 
systematic quantitative calculations of quarkonium processes, but it also leads 
to new qualitative insights into quarkonium physics. 
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